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Abstract 

k£ In this paper, we consider a problem for finding the locations of electromagnetic inhomo- 
geneities completely embedded in homogeneous two layered medium. For this purpose, we 
present a filter function operated at several frequencies and design an algorithm for finding 
the locations of such inhomogeneities. It is based on the fact that the collected Multi-Static 
Response (MSR) matrix can be modeled via a rigorous asymptotic expansion formula of 
the scattering amplitude due to the presence of such inhomogeneities. In order to show the 
effectiveness, we compare the proposed algorithm with traditional Multiple Signal Classifi- 
cation (MUSIC) algorithm and Kirchhoff migration. Various numerical results demonstrate 
that the proposed algorithm is robust with respect to random noise and yields more accurate 
location than the MUSIC algorithm and Kirchhoff migration. 
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1. Introduction 

In the non-destructive evaluation area, an inverse problem is finding the specific char- 
acteristics (location, geometry, internal constitution etc.) of small inhomogeneities from 
measurements of scattered or far-field data. This problem, which arises in fields such as 
physics, engineering, and biomedical science, is highly related to everyday human life and is 
still a challenging problem 0, @, 0] . 

Among them, finding the locations of anti-personnel mines (embedded in the soil) from 
electromagnetic data (measured in the air) is an interesting and important problem in the 
military services. The mines have different material properties from the surrounding medium 
and they are small relative to the area. The main purpose of this kind of application is to 
find their locations more accurately, not to retrieve complete information. For this purpose, 
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non-iterative MUSIC-type algorithm for finding locations of small inhomogeneities buried 
within a half-space at a fixed frequency has been developed, refer to jij. MUSIC type 
algorithms are advantageous in the sense that they are fast, stable, can easily be extended 
to the multiple inhomogeneities, and they do not require specific regularization terms that 
are highly dependent on the problem at hand. However, due to the reason that one is faced 
with an aspect-limited inverse problem in the reflection mode because sources and receivers 
are located upper half-space, sometimes inexact locations are identified. Hence, in order 
to obtain an accurate locations of such inhomogeneities, an alternative algorithm that can 
overcome these problems is necessary. 

Motivated by the above fact, we propose an effective, non-iterative location search al- 
gorithm at multiple frequencies that can work on limited view data in order to find the 
accurate locations of small electromagnetic inhomogeneities completely embedded within a 
homogeneous (lower) half-space. The starting point is that the collected MSR matrix can 
be approximated by a rigorous asymptotic expansion formula of the scattering amplitude in 
the presence of such inhomogeneities. This approximation leads us to proceed the singular 
value decomposition of the MSR matrix and to identify the structure of singular vectors. 
Applying the structure of singular vectors, a filter function operated at several frequencies 
can be designed for finding exact locations of inhomogeneities. Moreover, based on the 
statistical hypothesis testing, we can confirm that the proposed location search algorithm 
offers more exact information of locations than the established MUSIC-type algorithm and 
Kirchhoff migration. 

Various numerical examples with noisy data will illustrate the behavior of the proposed 
location search algorithm operated at several frequencies. Unlike the numerical setting in j^] , 
computational examples under the narrow/wide range of incident and observation directions 
configuration and the closely located small inclusions situation show the feasibilities and 
limitations of the proposed algorithm. 

This paper is organized as follows. In section [2J we briefly survey the two-dimensional 
direct scattering problem for two-layered medium and introduce an asymptotic expansion 
formula for the scattering amplitude. In section |3l we design a location search algorithm 
by producing a filter function operated at several frequencies and compare its detection 
performance with traditional MUSIC-type algorithm and Kirchhoff migration. In section 
HI corresponding numerical experiments with random noise are shown and compared with 
the result via MUSIC-type one and Kirchhoff migration to demonstrate its performance. In 
section [5] we give a brief conclusion. 

Finally, we will refer to 13J as a useful reference that investigates the so-called reciprocity 
gap MUSIC as another algorithm linked to the MUSIC one. In this reference, the author 
shows that reconstructed positions of the small inhomogeneities via reciprocity gap MUSIC 
algorithm are slightly more accurate than the traditional MUSIC one. 



2. Survey on two-dimensional direct scattering problem 

In this section, we briefly discuss the two-dimensional, time-harmonic electromagnetic 
scattering from a small inhomogeneity buried within a homogeneous half-space. A more 
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detailed description can be found in 

Let us decompose the two-dimensional space R 2 into the lower and the upper half-spaces 

as 

R 2 __ = {x= ( Xl ,x 2 ) T eM 2 : x 2 <0] and R^ = {x = (x x , x 2 ) T G R 2 : x 2 > 0} , 

respectively and assume that these spaces are homogeneous. Let <9R 2 denote the border 
between R 2 and R+, represented as 

<9R 2 = {{x u x 2 ) T G R 2 : x 2 = 0} . 

Throughout this paper, we assume that all the electromagnetic small inhomogeneities 
D m , m — 1, 2, • • • , M, are completely embedded in the lower half-space R 2 . For convenience, 
let T> denote the collection of such inhomogeneities as 

M M 

V= |J D rn = [J(z m + rB m ), 

m=l m=l 

where r is a small (with respect to the wavelength of the electromagnetic field in the em- 
bedding space at the given frequency of operation u) positive constant which denotes the 
diameter of the inhomogeneities, Bj is a simply connected smooth domain containing the 
origin, and Zj G R 2 , indicates the location of inhomogeneities. Throughout this paper, we 
assume that these inhomogeneities are separated enough from each other and from <9R 2 . 

All materials are characterized by their dielectric permittivity and magnetic permeability 
at the given frequency u; e_, e + and e m denotes the electric permittivity of R 2 , R^ and D m , 
respectively. The magnetic permeabilities //_, fi + and be defined analogously. Using 

these notations, we can define the piecewise constant electric permittivity < e(x) < +oo 
and magnetic permeability < /u(x) < +oo as 

( e + for xGl} ( fi+ for x G R^ 

e(x) = < e- for x G R 2 _\V and /x(x) = < p.- for x G R 2 _\V 
{ e m for x G D m [ p m for x G D m , 

respectively. For convenience, we also define the electric permittivity £o( x ) and magnetic 
permeability /io(x) when there are no inhomogeneities such as 

, , f e + for x G Ri. , , . f a + for x G Ri 

£o(x) = {s. for x G R 2 ! and MX) = \ a*, for x G R 2 ^, 

respectively. Accordingly, the piecewise positive real- valued wavenumber fc(x; u) reads as 



fc(x; cu) 



k + (u) = Uy/e + p + for x G R+ 
A;_(cj) = uj^/eZJIZ for x G R 2 . 



Let = (0i,8 2 ) T be a two-dimensional vector on the unit circle S l C R 2 and M; nc (x) = 
exp(ik + (u)6 ■ x) be a planar incident wavefield generated in the upper half-space on R 2 ^. At 



a given frequency u, «tot(x; co) denotes the time-harmonic electromagnetic total field which 
satisfies the following two-dimensional Helmholtz equation 



/i(x) 



VM to t(x;w) I + u 2 e(x)u tot (x;uj) = in M 2 , 



with transmission conditions holding at boundaries dM. 2 and dD m as 



[tttot(x; u)] 



1 du tot (x]UJ 



_//(x) 



on dM. 2 and [it tot (x; u;)] 



1 <9w to t(x;a; 



_/i(x) <9z/ m (x) 



on dD r 



respectively. Here, v m denotes the unit outward normal to dD m , m = 1,2, ■ ■ • ,m. 

Let Mbac( x ; a; ) be the solution to the Helmholtz equation ([!]) in the absence of inhomo- 
geneities. Then the scattering amplitude is defined as a function K : (S' 1 \{0}) x (S' 1 \{0}) x 
R — > C that satisfies 



{k-(u)} 2 n+(l+i) exp(ik(y;iu)\y\) 
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as |y| — > oo uniformly on 1? = A. 

In order to represent asymptotic expansion formula of K^fi, 0), we need some ingredients. 
First, based on the fact that the scattered field data acquisition is possible only on the upper 
half-space M 2 ,, we divide the unit circle S 1 into 



Si 



tins 1 and si 



ns 1 . 



Next, by letting £ = f±M define a vector 0(0; w) : £^{0} — ► C 2xl as 



T 



and a function $(0; u) : S' 1 \{0} — >■ C as 



//_f0 2 + // + sign(0 2 ) v /l-f 2 0? 



(2) 



(3) 



respectively. With these, an asymptotic expansion formula of scattering amplitude can be 
written as follows. This formula plays a key role of the location search algorithm that will 
be designed in the next section. 

Theorem 2.1. For every •& 6 S\ and G Si, the asymptotic formula for the scattering 
amplitude K(i9, 0) at frequency u is expressed as 



M 



m=l 



K{-&, 0; u) =r 2 $(tf ; u)${0; u) £ l E , m \B m \ + J^^fr, u) ■ P(z m ) • 0(0; u) x 



exp ( - ik_(uj)[4>(&,u) - 0(0; u)] ■ z m ) + o(r 2 ), 
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(4) 



where the remaining term o(r 2 ) is independent of i? G S], 6 G Si and the set of points 
{z m } m=1 , P(z TO ) is a 2 x 2 positive, symmetric matrix 



At- + /Um 

and constants ^^^m and ^ e ,m are given by 

le,m = 1 and 7 M , m = I. 

3. Non-iterative location search algorithm &c its performance 

3.1. Non-iterative location search algorithm at multiple frequencies 

We apply the asymptotic formula for the scattering amplitude (|1J) in order to build up 
the location search algorithm. For this purpose, we will use the eigenvalue structure of the 
Multi-Static Response (MSR) matrix K := K(u) = (Kji(u)) G C N ° h * xN ™, whose element 
Kji{u) := K(-dj, Of, io) is the scattering amplitude collected at observation number j for the 
incident wave numbered I. In this paper, we denote 

{fy : j = 1, 2, • • ■ , N obs } and {0, : I = 1, 2, ■ ■ ■ , N inc } 

be the set of observation and incident directions, respectively. Note that when the upper half- 
space is more refractive than the lower one, i.e., if k+ > k-, the number N oha of propagating 
transmitted waves might be less than N inc (see [3] for instance). Then since jl— th element 
of the MSR matrix Kji(u) can be approximated as 

Kji{u}) =K('d j ,O l ;co) 

M , 

^r 2 $(tf j; u)§{Ou u) he, m \B m \ + 7^0(0; u) ■ P(z m ) ■ 0(0 

7Tt=l ^ 

exp f-ik-(uj)[(l)('& j ;uj) -0(0j;w)] • z m Y 



(5) 



MSR matrix K can be decomposed as follows: 

X = OEHF, (6) 
where E G C 3Mx3M is a diagonal matrix with components 



E 



E e 
E 



for 



E e = M x M diagonal matrix with components r 2 7 £jm |.B m |, 
E M = 2M x 2M diagonal matrix with 2x2 blocks - r 2 7 Mim P(z„ 
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matrix D G C^^^ is of the form 
with iV obs x 1 matrices 



5f 



u) exp ( — ik^(u)4>('& 1 ; u) ■ z r 



m I 5 j 



j2(to-1)+s 



e s • w) exp ( — ik_(u)4>('&i; u) ■ z 



ml; ) 



e s ■ 0(^JV ob8 ; w)$(0jv obs ; w) exp - ifc_(w)0(0w ob .; w) ■ z 



and matrix EI e c^incxSM can ^ e wr j^ en as follows 
with N inc x 1 matrices 



Hf Hi e 2 



HI' 



jj2(m-l)+s 



w) exp I ik_(tu)(p(6 1 ; u) • z r 



ml) j 



$(0jv inc ; w) exp ik_(uj)(j)(0 N . nc ; u) ■ z r 



e s ■ <j>{Oi\ oj)Q{di] u) exp I ik^(u)4>(0i, to) ■ z 



ml) ) 



e s • 4>{0N inc ; u)®(0 N . nc ; exp ( ik_(uj)(j>('d Ninc ; u) ■ z 



Here {e s } s=12 = {ei = (l,0) T ,e 2 = (0, 1) T } is an orthonormal basis of 1R 2 . Based on the 
decomposition ([6]), a location search algorithm can be established as follows. 

1. (Singular Value Decomposition) Let us perform Singular Value Decomposition (SVD) 
of matrix K and let M be the number of nonzero singular values for the given u. Then, 
K can be represented as follows: 



K = U(w)S(w)V*(w) 



M 

E 

m=l 



where superscript * denotes the complex conjugate, cr m (w) are the singular values, 
U m (u) and V m (w) are the left and right singular vectors of IK, respectively for m = 



1,2,.-. ,M. 
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2. (Structure of singular vectors) For test vectors Cd, Ch € M 3 \ {0}, define vectors 



d : 



^AT obs Xl an( J ; 



X 



£_ N inc X 1 



as 



dfx; u) 



c d ■ (f , 0(#i; w) exp ( — ifc_(o;)^)(i? a ; w) ■ x 



c d • (1, 0(#iv obs ; w))$(^7v obs ; w) exp - ik-(u)4>(0 NobE ; w) • x 



(7) 



and 



hfx; a;) 



c h ■ (1, 0(#i; o;))$(0i; w) exp ( ik-(uj)<f>(0i; uj) ■ x 



Ch • (1, 0(^jv inc ; w))^(^JV inc ; w) exp Uk-(u)<l>(0 Niae ; uj) ■ x 
respectively. With this, generate corresponding normalized unit vectors 

W ( \ d ( X ^) A W I \ h fx;U;) 

W D (x; oj) = — — and W H (x;u;) 



1 t 



dfx; uj) 



hfx; 



Note that the structure of vectors Wd(x; uj) and Wh(x; uj) is motivated by the matrix 
D and HI in (jBj), respectively. Then by virtue of [14(, following relationship holds for 



m = 1,2, ••• ,M, 



U m (u>) ~ W D (z m ;w) and V m (w) ~ W H (z m ; w). 



(9) 



Since the first M columns of the matrix U(u) and V(w), {Ui(u;), U 2 (o;), ■ ■ • , XJ m (uj)} 
and {Vi(u), V 2 (a;), ■ • • , V m (u;)}, are orthonormal, we can observe that 



W D (x;u;)*U m (u;) ^ and W H (x; u)*V m (u) ^ 0, if x = z m 
W D (x;u;)*U m (u;) « and W H (x; w)*V m (w) « 0, if x ^ z m , 



(10) 



for m = 1,2, • • • , M. 

3. (Filter function) For a search domain Q C Ml, construct a normalized filter function 
F:flxN — > M at several frequencies {w/ : / = 1, 2, • • • , F} as 



F(x; F) = — 



^(w D (x; U/ )'U m J (w H (x; 07)^(0;,)) 

f=l m=l ^ ' ^ ' 



Then, based on the observation (flOj) . we can find locations z m G -D m by finding x 
which satisfies F(x; F) ml. 

Remark 3.1 (A priori information). In order to build up the filter / T77]) . we need a priori 
information of the values e_ and //_. If one has no information of or simultaneous 
reconstruction of multiple parameters must be performed (see ll l\ . Section 10.3] for instance). 



The location search algorithm is summarized in Algorithm [TJ 
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Algorithm 1 Location Search Algorithm 



procedure LSA(F) 

identify values e_ and /x_ 
initialize K(pc) 
for / = 1 to F do 

collect MSR matrix data K(u;/) 
perform SVD of K(tuj) 
discriminate nonzero singular values 
choose U m (uf) and V m (o>/) 
for x e f] C 1 2 _ do 

generate Wd(x; w/) and Wh(x; Uf) 
initialize 7(x, /) 
for m = 1 to M do 

J(x,/) <- J(x,/) + (W D (x; Uj )«U m ( U/ ))(W H (x; W/ )*V m ( Wj )). 
end for 
A(x)^7(x,/) 
end for 
end for 

plot F(x;F) = \K(x)\/F 



find x = z m G D r< 
end procedure 



3.2. Some properties of normalized filter function F 

At this moment, we explore some properties of normalized filter function F in (|TT|) . For 
this purpose, we assume that there exists one inhomogeneity (m = 1) and 

c d -(l,<t>{O n ;u f )) ^0, Cil -{l,<t>{O n ;u; f )) ^0 and $(0 n ; W/ )^O (12) 

for all n = 1, 2, • • • , N and / = 1, 2, • • • , F. Then applying relation (j^J) to (JTT1) yields 



F(x; F) = — 



(wDixiw^Uj^j IwHixiu/j^^J 



-|F D (x;F)F H (x;F)| 



where 

F D (x;F) = J] 

n=l 
A 7 " 

F H (x;F) = ^ 



A? 



71=1 



|c d - (lX^;^))^^;^)!- 
|d(z m ;a; / )||d(x; 

|c h - (l,0(6> n ;u; / ))<l»(6> w ;a; / )| 2 
|h(z m ;w / )||h(x;u; / )| 



exp ^ - ik-(uf)^(-& n ;u}f) ■ (z m - x) 
exp I ik-(uf)<j>(d n ; u f ) ■ (z m - x) ) . 
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With this, by letting # n := (i?" , fi^Y \ '■= (#1, #2 ) T an d z := ( z i, z 2) T , we can observe 
that 

JV 

F- 



D (x; F) ~ f>xp ( - zfc_( W/ ) (^(^i - *i) + V 7 ! - (W(^2 - x 2 )) ] 

n=l ^ ' 

F H (x;F) ~ ^Texp (^^(^ - x x ) - ^Jl - Wi 

n=l ^ 



Z 2 ~ X 2 , 



where A ~ means that there exists a constant C such that A = BC. Throughout 
this paper, we assume that iV obs and iV inc are even number iV obs = 2L, {fl n } and {#„,} are 
symmetric to y-axis, i.e., if tf n := then = _tf2£-n-i and = ^p-n-i^ and go 

on. Then we can explore some properties of F(x; F) as follows 

1. Assume that k+(ojf) satisfies k+(cj/) > k-{uf) and£$" > 1. Then since y/l — (£$i) 2 G 
C\{0}, by letting a/1 — (^") 2 = «Pn, Fd(x; F) can be written as 

F D (x; F) ~ ^exp ^ - iA; + (u;/)£t?i(zi - zi)^ exp ^p n /c_(a; / )(2:2 - s 2 )^ . 
This means that F(x; F) ~ 1 when 

S7T S7T 7T 

Xr = Zr + - — : — -— , Xt — Zi + - — ; — — — I — and x 2 = z 2 
k+(cj f )6f 2 

for all s G Z, / = 1,2, •■■ ,F, and n = 1,2, *•• ,N obs . This relation tells us that 
in order to obtain an accurate location z, one must applies high frequency and large 
number iV obs . Nevertheless, some replicas will appear along the x— axis. However, if 
one adopt symmetric observation configuration, Fd(x; F) becomes 



F D (x; F) ~ cos ( - ik +("f)W 

n=l ^ 



(zi - Xi)j exp \p n k-(uf)(z 2 - x 2 )J . 

Hence, we can obtain more accurate location z. See Figure HI 

2. When the value F is large enough, due to the different values of Uf, F(x; F) will yields 
more accurate location z. This means that application of multiple frequencies will 
enhance detection performance. 

3. Assume that k + (uf) ~ k-(ujf) and i?" — >■ ±1 for some n, i.e., one has wide observation 
direction, then since yl — (<f$i) 2 ~ 0, location x such that F(x; F) ~ 1 is independent 
to the z 2 . Hence we cannot identify location z via F(x; F), refer to Figure [TU1 

4. If fc+(w/) < k-(u f ) then since ^1 - (£^) 2 G K\{0}, F(x; F) « 1 when 

M"/) " zi) " V 7 ! " (^?) 2 (^2 - ar 2 )) = svr or svr + | (13) 

for all s G Z, / = 1, 2, • • • , F, and n = 1, 2, ■ • • ,N. Therefore, map of F(x; F) will 

offers exact location z but unexpected some ghost replicas will obstruct it. 
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5. If k+(uf) <C k-(uf), (|T3|) becomes 

fe_( W/ ) (e^(^i - an) - Vl - W(^2 - x 2 )) 



k+( Uf W(zi - an) - yJk-(Uf) 2 -(k + (Uf)^) 2 (z 2 - x 2 ) 

7T 

fc + (u>/)i?i(;si - an) - k-(cjf)(z 2 - x 2 ) = stt or S7r + - 



for all s e Z, / = 1, 2, • • • , F, and n = 1, 2, • • • , AT. This shows that F(x; F) » 1 if 

S7T (2s + 1)7T 

Si = 2i, x 2 « z 2 + fc ^ ^ or xi = z x , x 2 « ^ 2 + 2fc ^ ^ ■ 

Since k-(uf) is large enough, huge amounts of replicas will appear along the axis 
in the map of F(x; F) (see Figure [9]). 

The case of Fd(x; F) can be handled in similar manner. 

3.3. Introduction to MUSIC and Kirchhoff migration 

By combining the results of [H, 20], we can design a MUSIC-type image function at a 
single frequency u. We define a projection onto the null (or noise) subspace P no ise as 

P noisc (d(x; u)) = U m MU;»d(x; u>). (14) 

m>M 

Then the image of x = z m , m = 1, 2, • • • , M, follows from the computation via a MUSIC- 
type imaging functional Fmusic : I^ 2 - x K — > 

F MU sic(x;a;) = — — — . (15) 

|P noise (d(x;u;))| 

With this, we can find locations x = z m which satisfy F MU sic( x ; u>) = oo. 
We introduce the traditional Kirchhoff migration Fkir : M 2 x R — > K 



F K ir(x;o;) = |W D (x; U )W H (x;o;) 



(w D (x; w)*U m (a;)) (w H (x; u;)*V m (^)) 

(16) 



Then similar to the filter function (ITTj) . map of Fkir(x; cj) will yields locations z m . 

3.^. Comparison of detection performance 

Now, we will briefly compare the detection performance of (|15p and ffTB]) . Roughly 
speaking, the following relationship holds for F > 1 

Fmusic (x; w) < F K ir(x; w) < F(x; 1) < F(x; F), (17) 

where A < B means B offers more accurate location z m G D m than A. 
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First, based on the recent work |3|, the relationship Fmusic( x ; <j Fkir(x; uj) holds for 
homogeneous space case and we can easily verify that it also holds for two-layered medium 



problem interested herein (see [18J for numerical experiments). 

Next, when K is affected by random noise, significant changes of singular values will 
appear. In this problem, such noise generates many nonzero singular values (see 0, Figure 
5.25] so that (1T6"|) generates a result with poor resolution. However, ( ITT]) is not influenced 
by the singular values and therefore, map of F(x; 1) yields a better result than Fkir(x;o;), 
refer to Figure El 

Remark 3.2 (Synthetic Aperture Radar (SAR) and Kirchhoff migration). Synthetic 
Aperture Radar (SAR) is one of the classical back-projection imaging technique developed for 
airborne radar applications (see |J& JL5,[llJ/ and references therein). Based on the research 



in 



'23] , results via SAR and Kirchhoff migration are almost the same. 



Finally, based on the section l3~2l we can examine the relationship F(x; i r i)<F(x; F?) holds 
for F\ < Fi. It is worth emp hasizing that we can observe this relationship via statistical 
hypothesis testing (see 

IBIS for detailed discussion). 



4. Numerical simulations and discussion 

In this section, various numerical results are presented to demonstrate the effectiveness 
of the proposed algorithm. Same as the numerical configuration in [H, we choose three small 
homogeneous inhomogeneities, D\, D2 and D%, embedded in the lower half-space. They are 
taken as ball of radius r = 0.1 and are centered at zi = (0.63, —2.47), z 2 = (1.72, —4.97) and 
z 3 = (—2, —3.63), respectively. The applied frequency is Uf = jj, where A/, / = 1, 2, • • ■ ,10, 
are given wavelengths. In this paper, frequencies Uf are equi-distributed within the interval 
[2tt, 1^]. The observation and incident directions i9j and 0\ are taken as 

7T (j — 1)7T 7T (I — 1)7T 

Vj = (cos Q, sin Q) , Q = T+TTTTT 7T> an d &i = ~ (cos sin ^) , q = -+ 



4 2(iV obs - 1) ' ^ ^ J ' ^ 4 2(N inc - 1) 

respectively for j = 1, 2, • • ■ , iV obs and I = 1, 2, ■ ■ • , N inc . See Figure [TJ for an illustration of 
the test configuration. 

Throughout this section, we adopt the squared search domain Q = [—3, 3] x [—6, 0] C R?_. 
The step size of the search points x 6 Q is taken of the order of 0.05 and vector c of (j7|) is 
selected as 

1. Permittivity contrast case: c = (1,0,0), 

2. Permeability contrast case: c = (0,1,5), 

3. Both permittivity and permittivity contrast case: c = (5, 1, 1). 

Note that this selection of c satisfies (I12p. A detailed discussion about the choice of the vector 



c can be found in 2l|, Section 4.2.1]. In every example, the data set of the MSR matrix is 

computed within the framework of the Foldy-Lax equation, refer to I24J 1. Then, a 

white Gaussian noise with 20dB SNR(Signal to Noise Ratio) is added to the unperturbed 

11 



Figure 1: Sketch of the test configuration. 



data in order to show the robustness of the proposed algorithm. Note that from the various 
numerical experiments in 0, 18, 20 1, similar results were obtained for both cases e + > e_ 



and e + = e_ (permeability and both contrast cases too). Thus, we do not consider the cases 
e + = E- and/or //+ = //_. 

4-1. Permittivity contrast case: e m ^ and \i m = //_ = //+ 

At this stage, we consider the purely dielectric contrast case. In this case, we set /x(x) = 1 
for x G M 2 . For e + > case, we choose the values e + = 5 and = 4 and permittivities 
e m of D m equal to 2, 5, 3 for m = 1,2, 3. As already mentioned, the number of observation 
directions N ohs must be smaller than the number of incidence directions N inc . Hence, we 
choose N inc = 10 and iV obs = 6 directions. 

Let us check the detection performance of ffTTj) and ffTET) . As the results illustrated in 
Figure [U map of F(x; 1) offers more accurate location of small inhomogeneities (specially, 
D3) than the one of Fkir(x;o;). Moreover, by comparing the result in [J, FIG 5.3] and the 
map of Fkir(x; u>), we can see that the relationship Fmusic( x ; u>) < Fkir(x; oj) holds. 

Figure [3] shows the influence of the number of applied frequencies. By comparing maps 
of F(x;F), we can easily observe that large number of F guarantees an exact location 
of D m . In various numerical tests, we observed that if one applied more than F = 7 
different frequencies, map of F(x; F) yields an accurate location of D m , so we adopt F — 10 
frequencies in order to guarantee an admissible result. 

In order to examine the influence of direction configuration, odd number of incident and 
observation directions is applied and corresponding results are illustrated in Figure HI By 
comparing result in Figure [3J as we mentioned in section 13.2} some replicas appeared under 
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Figure 2: (Compare the detection performance) Map of Fkir(x) (left) and F(x; 1) (right) for uj = 2tt when 
e+ = 5 and e_ = 4. 



the odd number configuration but when the number of directions and frequencies is large 
enough, we can identify location of D m accurately. 

For e + < E- case, we choose the values e + — 1 and 6- = 5 and permittivities e m of 
D m equal to 2, 4, 3 for m = 1, 2, 3. In this case, although a few ghost replicas appeared, the 
location of three inhomogeneities are successfully identified, refer to the right-hand side of 
Figure EJ 

4-2. Permeability contrast case: e m = e_ = e + and /i m ^ n_ 

Now, let us consider the purely magnetic permeability contrast case. In this case, we set 
e(x) = 1 for x G M 2 . Similar to the section I4.1[ we choose the values fi + = 5 and //_ = 4 
and permeabilities \i m of D m equal to 2, 5, 3 for m = 1, 2, 3. We apply F = 10 frequencies, 
and iV inc = 14 and iV obs = 8 as directions of incidence and observation, respectively. For 
fi + < H- case, we choose the values /i+ = 1 and = 5 and permeabilities fi m of D m equal 
to 2, 4, 3 for m = 1, 2, 3. Figure [6] shows the corresponding result. Similarly to Figure the 
location of three inhomogeneities are successfully identified. 

4-3. Both permittivity and permeability contrast case: e m ^ s_ and \i m ^ /i_ 

In this case, we consider both permittivity and permeability contrast case. Three different 
situations of interest are considered: 

• e + 7^ E- and = 

• and /i + 7^ 

• e + 7^ e_ and /x + 7^ 
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Figure 3: (Influence of number of applied frequency) Map of F(x; F) for F = 1 (top, left), F — 3 (top, 
right), F = 5 (bottom, left), and F = 7 (bottom, right) when e + = 5 and e_ = 4. 



First, let us consider the case e + > £_ and /i + = //_. We set /x(x) = 1 for x G M 2 , e + = 5 
and e_ = 4 and permittivities e m and permeabilities fi m of _D m equal to 2, 5, 3 for m = 1, 2, 3. 
We apply F — 10 frequencies and adopt iV inc = 20 and iV obs = 12 as directions of incidence 
and observation, respectively. The result, as illustrated in Figure UJ remains very good and 
should be acceptable. Now, let us consider the case e + < e_ and /x + = //_. Let e + = 1 and 
£_ = 5 and permittivities e m and permeabilities /i m of D m are equal to 2, 4, 3 for j = 1, 2, 3 
while keeping remaining test configurations. Similarly to the previous example, a good result 
appeared, refer to Figure [71 With a similar argument, we can obtain a good result when 
e + = 6- and ji + ^ //_, refer to Figure [8j 

For the case e + ^ e_ and //+ 7^ we set e + = /x + = 5, e_ = /i_ = 4 and £ + = ji + = 1, 
= /i_ = 5 while keeping the configuration of the previous situation. The result is 
exhibited in Figure |9j Although we could obtain a reasonably good result when e + > gr- 
and fi + > one can not determine accurate locations of inhomogeneities when e + < gr- 
and /i + < based on the discussions in section 13.21 
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Figure 4: (Influence of direction configuration) Map of F(x; F) for N obs = 5, N- inc — 7, F = 5 (left), and 
N obs = 11, iV inc = 13, F = 10 (right) when e + = 5 and e_ = 4. 




Various results in this paper show that the proposed algorithm is very stable and effective 
but it still contains some factors for further improvements. 

4-4- Robustness with respect to random noise and influence of range of incident and obser- 
vation direction 

At this moment, we add a white Gaussian noise with lOdB SNR to the unperturbed data 
and change the range of incident and observation directions. For this purpose, we adopt the 
same test configuration as section 14.31 except that the ranges of observation directions 



&j = ~ ( cos 0, sin Q) , j = 1, 2, ■ ■ • , N obs , 
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Figure 6: (Permeability contrast case) Map of F(x; 10) when fx + 
(right) 



5, fi- = 4 (left) and /i+ — 1, /x_ — 5 
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Figure 7: (Both permittivity and permeability contrast case) Map of F(x; 10) when e + = 5, e_ = 4, //(x) = 1 
(left) and e + =!,£_ = 5, //(x) = 1 (right). 



are varied by changing the values of (i and 0v O bs- Incident directions are also changed 
similarly. Figure [10] shows the map of F(x; 10) when e + = /i + = 5, e_ = /i_ = 4 under the 
narrow and wide range of observation and incident direction configurations. Regarding the 
top, left-hand side of Figure [TUJ we can see that under the narrow range of incident and 
observation direction configuration, map of F(x; 10) contains some ghost replicas similar to 
the case of e + < e_ (see Figure [7]). From the top, right-hand side of Figure [101 we can clearly 
identify three locations of z m . Therefore, we can conclude that the proposed algorithm is 
robust with respect to the large amount of random noise since the normalized filter function 
( II ip is not significantly influenced by the noise. Unfortunately, based on the section I3.2[ by 
regarding the bottom line of Figure [TU] map of F(x; F) yields poor result when the range 
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Figure 8: (Both permittivity and permeability contrast case) Map of F(x; 10) when e(x) = 1, /i + = 5, 
/i_ = 4 (left) and e(x) = 1, /x + = 1, /i_ = 5 (right). 




of directions become wider. Hence, we can conclude that when the range of incident and 
observation directions become wider, we cannot find the location z m via map of F(x;F). 
This shows a limitation of the proposed algorithm. 

4-5. On the Rayleigh resolution limit 

Now, we briefly consider the image resolution. From the Rayleigh resolution limit, we 
can distinguish any inhomogeneities D and D' when 

dist (£>,£>') > ~. 
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In order to examine such phenomenon, we consider the detected location of two disks D± 
and -D5 of the same radius r = 0.01 with permittivity contrast case. The centers of -D4 and 
D 5 are selected as z 4 = (—0.1,-2.5) and z 5 = (0.1,-2.5), respectively. In Figure HH we 
illustrate corresponding results with large and small wavelengths while keeping remaining 
test configurations. Therefore, it is hard to distinguish two inhomogeneities from the image 
via map of F(x; 10) with large wavelengths. On the other hand, when we apply smaller 
wavelengths (i.e., higher frequency) than the previous one, the two inhomogeneities become 
distinguishable. 

5. Concluding remarks 

In this paper, we suggest a location search algorithm operated at several time-harmonic 
frequencies in order to find accurate locations of small electromagnetic inhomogeneities 
completely embedded within a homogeneous lower half-space. The approach is based on 
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Figure 11: (Permittivity contrast case with e+ = 5 and e_ = 4) Map of F(x; 10) when Ai = 2, Aio = 1 (left) 
and Ai = 0.2, Am = 0.1 (right). 



the asymptotic formulation due to the existence of small electromagnetic inhomogeneities. 
Throughout numerical simulations, we can conclude that the proposed algorithm not only 
performs quite well even in the existence of random noise but also successfully improves 
existing limitations (poor longitudinal resolution against an excellent transverse resolution) 
of MUSIC algorithm proposed in [i| and Kirchhoff migration. In addition, it still has 
some points of improvement for finding the locations of inhomogeneities under the situation 
e + < e_, n + < ji- and wide range of incident/observation directions. 

It is worth mentioning that such results obtained at low computational costs can be a 
good initial guess of a level-set evolution 0,0, H, 22[ or of any other standard iterative 
algorithm. Although only two-dimensional problem have been considered herein, we expect 
that the proposed strategy, e.g., asymptotic formula, filter design, etc., could be extended to 
the three-dimensional problem, refer to 0, [l5| for related works. Moreover, inconveniences 
of the proposed algorithm are also the same as those of SAR0- Hence, comparison of SAR 
and the proposed algorithm will be an interesting subject. 

Finally, we would like to emphasize that the proposed algorithm can be extended to the 
shape identification of electromagnetically thin, arc-like, penetrable inhomogeneities, refer 



to 



171 . Il8l . 120] . Although, further mathematical investigation is necessary, we believe that it 



can also be extended to the identification of small or extended perfectly conducting cracks. 
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